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Abstract 

The recently introduced self-organizing height-arrow (HA) model is numeri- 
cally investigated on the square lattice and analytically on the Bethe lattice. 
The concentration of occupied sites and critical exponents of distributions of 
avalanches are evaluated for two slightly different versions of the model. The 
obtained exponents for distributions of avalanches by mass, area, duration 
and appropriate fractal dimensions are close to those for the BTW model, 
which suggests that the HA model belongs to the same universality class. For 
comparison, the concentration of occupied sites in the HA model is calculated 
exactly on the Bethe lattice of coordination number g = 4 as well. 
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I. INTRODUCTION. 



The study of different cellular automata, which exhibit Self-Organized Criticality 
(SOC) m, has been a subject of great interest in recent years. These models serve as 
tractable limits of real dynamic systems with many spatial degrees of freedom, in which one 
might hope to gain understanding of possible mechanisms of SOC. Unfortunately, in most 
cases our current knowledge of the effects of SOC is rather limited. 

The main peculiarity of the dynamic dissipative models which are self-driven into the 
SOC state is the presence of the power-law in distributions of quantities such as avalanche 
mass, duration, etc. These distributions are characterized by a set of exponents. One of 
the most intriguing questions concerns the classes of universality of these models. There 
are several attempts to shed light on this problem Recently, Ben-Hur and Biham 

proposed a classification scheme for the 2d models both stochastic and deterministic. They 
found that the original Bak, Tang, Wiesenfeld (BTW) model belongs to the universality 
class of undirected models, directed models form a separate class, and the two-state Manna 
model 0] belongs to the universality class of random relaxation models. Later on, Nakanishi 
and Sneppen examined several Id sandpile models and suggested that the two-state 
Manna model and rice pile model belong to the same universality class. 

This classification scheme is based on the type of relaxations at each site of the lattice. 
In the BTW model the particles from the toppled sites are uniformly redistributed among 
its nearest neighbors, whereas in the two-state Manna model the set of neighbors is chosen 
randomly. It is possible to introduce more complicated dynamical rules of relaxations at each 
site of the lattice that will depend on some period T where T is the number of topplings. 
During this period after each toppling the redistributions of particles from the given site 
form a minimal nonperiodic sequence. The BTW model can be considered with the period 
T = 1. In the two-state Mannna model the sequence of topplings at each site is stochastic 
without any periodicity, therefore, one can put T = oo for this model. 

In this paper, we investigate the recently introduced self-organizing height-arrow (HA) 
model 0J^. It combines features of the BTW model 0,10] and self-organizing Eulerian 



walkers model (EWM) ||7|JTT[|. The model is a cellular automaton defined on any connected 
undirected graph. In this model, each site of the graph can be occupied by a particle or 
can be empty. Addition of the particle to the occupied site makes it unstable and causes its 
toppling. The site becomes empty and the particles are transferred to the nearest neighbors. 
The redistribution of particles from an unstable site is governed by the second site variable, 
an arrow. Each outgoing particle from the toppled site turns the arrow by the prescribed 
rule and the new direction of the arrow determines the destination point for this particle. 

In the HA model T is formed by the nonperiodic sequence of turns of the arrow at each 
site of the lattice. For simplicity, it is convenient to assign the same period T for all sites 
of the lattice. Thus, one can define the HA model with increasingly complicated dynamical 
rules. These pseudo-random models tend to the random ones for large T — > oo. 

The goal of the paper is the study of the HA model with T = 2, when two topplings 
restore initial direction of the arrow. The model evolves at long times into the steady state 
which is identified with the SOC state as distributions of dynamic characteristics of the 
model show a power-law behavior. The obtained exponents for two slightly different types 
of this model are very close to those for the BTW one, which suggests that the HA model is 
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in the universality class of undirected models. We also study the main static characteristic of 
the model, the time averaged density of occupied sites. This quantity is obtained numerically 
on the square lattice and exactly on the Bethe lattice of coordination number g = 4. 

The paper is organized as follows. The definition of the model with two slightly different 
types of evolution rules is presented in the next section. Sec. ^ is devoted to the numerical 
investigation of the HA model on the square lattice. Extrapolating results from finite size 
lattices we estimate the concentration of occupied sites in the model. The critical power- 
law exponents and scaling relations among them are defined in the framework of finite-size 
scaling analysis. We present results for the values of these exponents in the SOC state for 
distributions of various quantities of the HA model. Then, in Sec. |IV|, we exactly calculate 
the concentration of occupied sites on the Bethe lattice of coordination number q = 4. 
Discussion and conclusions are presented in Sec. |V|. 

II. THE SELF-ORGANIZING HEIGHT-ARROW MODEL ON THE SQUARE 

LATTICE. 

The HA model we are going to investigate is defined as follows. To each site i of the 
two-dimensional Lx L square lattice is assigned a height variable Zi G {0, 1, ...} and an arrow 
directed north, east, south or west from i. We start with an arbitrary initial configuration of 
heights and arrows on the lattice. Initially, we drop a particle on the randomly chosen site i. 
The succeeding evolution of the system is determined by the following rules. We increase the 
height variable at the site i by 1, Zi ^ Zi + 1. If the site i is already occupied by the particle, 
it topples {zi Zi — 2). To redistribute the particles from the toppled site i among its 
nearest neighbors, we turn the arrow twice according to the prescribed rules. For the given 
period T = 2, there are only two non-equivalent sequences of turns of the arrow at the given 
site which preserve the model from being directed. Hereafter, these sequences of turns will 
be distinguished as N-E-S-W-N and N-S-W-E-N. After each turn the new direction of the 
arrow points to the neighbor sites to which we will transfer particles at the next time step. 
This process continues until a stable configuration is reached. The sequence of topplings of 
unstable sites forms an avalanche which propagates through the lattice. After an avalanche 
ceases, we go on by adding a new particle and so on. 

A given configuration of the model is a set of directions of arrows and heights. The total 
number of them is 8^^^. During the evolution of the system the arrow at any site might 
be only in two positions due to the fact that the two subsequent tophngs of the site restore 
the initial position of the arrow. Therefore, the set of configurations of the model falls into 
2-LxL equivalent classes which are determined by initial configurations of arrows. 

Starting from a certain configuration of arrows and an arbitrary configuration of occupied 
sites, the model evolves through transient states into a dynamic attractor which is critical. 
This attractor is identified with the SOC state as different dynamic characteristics of the 
model show power-law tails in their distributions. The model being in the SOC state passes 
from one allowed stable configuration to another by avalanche dynamics. This critical state 
has been investigated in detail by Priezzhev . He defined operators corresponding to addi- 
tion of a particle at a randomly chosen site and showed that they commute with each other. 
The algebra of these operators is used to calculate the number of allowed configurations of 
a given class in the SOC state. This number is shown to be equal to the determinant of the 
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discrete Laplacian matrix A of the square lattice. To check the given configuration to be 
allowed in the SOC state the modification of the burning algorithm was also introduced [§. 

III. NUMERICAL RESULTS. 

In order to investigate the static properties and avalanche dynamics in the HA model, 
we have made numerical simulations with high statistics. We consider square lattices of size 
L X L with open boundary condition and L ranging from 100 up to 600. The HA model has 
been studied for two different types of dynamics (N-E-S-W-N and N-S-W-E-N) of turns of 
arrows and various initial conditions. 

Starting from an arbitrary distribution of occupied sites and certain initial directions of 
arrows, the finite system evolves into a stationary state. In this state we have measured the 
time averaged density {p{z = 1)) and critical exponents for distributions of avalanches by 
mass (s), area (a) and duration (t). The mass s is defined as a total number of topplings 
in an avalanche whereas the area a is defined as a number of distinct sites visited by an 
avalanche. The simultaneous topplings of different sites in an avalanche at a given time is 
considered a single time step. The duration t is the number of this type of steps. For a 
more detailed description of the structure of an avalanche it is also useful to define a linear 
extent (diameter) of the avalanche cluster via a radius of gyration (r). We also measured 
the corresponding fractal dimensions 'jxy, where x,y = {s, a, t, r}. 

As is shown in Fig. |1|, the steady state is reached by the model after about 100 000 
avalanches on the square lattice of the linear size L = 600 for the system which is initially 
empty and with a random initial distribution of arrows. In this simulation, we were recording 
the averaged density {p{z = 1)) of occupied sites at each time when the avalanche ceases. 

As has been mentioned in Sec. |n|, the number of configurations of the model falls into 
2^^^ classes depending on the initial configurations of arrows. In our simulations of the 
HA model with the N-E-S-W-N dynamics we started from random initial configurations of 
arrows. Whereas, for the N-S-W-E-N dynamics the arrows were initially directed only east 
or south. The later case was chosen to simulate the scattering of particles at each toppling 
by 180° angle. 

Fig. 1^ displays the results of simulations for the time averaged density {p{z = 1)) in the 
stationary state. They slightly depend on the lattice size L and are well described by the 
equation {p{z = = Pc + cL~^. The numerical extrapolation of the L ^ oo limit gives 
the values for the averaged density: pc = \imL^oo{p{z = = 0.721 ± 0.001 (N-E-S-W-N 
dynamics) and Pc = 0.755 ± 0.001 (N-S-W-E-N dynamics). These values are a little higher 
in comparison with the stochastic two-state Manna model (see Table |ip. 

The form of avalanches in the HA model has a layered structure. A typical avalanche is 
shown in Fig. ^ where the number of relaxations in each site is marked by different scales 
of gray color. The sites with the same number of relaxations form a layer or shell. We 
have observed that layers group in pairs. In each pair a larger layer is a connected cluster 
with holes whereas a smaller one is a disconnected cluster without holes. Therefore, in the 
avalanche cluster there are a very few holes only near the surface in the first layer where 
each site topples once. 

In Figs. |, we present the directly measured distributions of avalanches by mass s, area a 
and duration t in a double logarithmic plot for the N-E-S-W-N dynamics and the lattice of 
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size L = 600. These distributions display a power-law behavior up to a certain cutoff which 
depends on the system size L. Since our simulations are limited by the lattices of finite size 
we ought to apply the finite-size scaling analysis |I2| . [T3| assuming the distribution functions 
scale with the lattice size L 



P{x,L) = L~^-Ux-L- 



(3.1) 



where fx{xL~'^'^) is a universal scaling function, x stands for s, a, t or r, and (3x and Ux 
are critical exponents which describe the scaling of the distribution function. The finite-size 
scaling ansatz ( |3.1| ) can be rewritten in the following form |T^: 



P{x,L) 



fx{x ■ L' 



(3.2) 



Let us suppose that distribution functions in the thermodynamic limit (L oo) show 
pure power-law behavior for large enough stochastic variables (s,a,t,r) 

P{x) ~ x"""^ , X > 1 , (3.3) 

where x G {s,a,t,r} are critical exponents. This conjecture is mainly supported by 



computer simulations and different heuristic arguments Therefore, comparing ( p.2| 



and ( |3.3| ) we get the scaling relations among these exponents 

rx = — ■ 

From the fact that (s) ~ in the undirected BTW-type models 
tional scaling relation [14 



(3.4) 

one can get an addi- 



Ps{2 - r,) = 2 . (3.5) 
If we also assume that the stochastic variables s, a, t, r scale against each other, the appro- 



priate fractal dimensions '~^xy can be defined via the following relations [15 



s ~ a 



7s a 



s ~ t^" 



7ar 



a ^ r 
t ~ r"^*-- 



(3.6) 



where '-^xy = lyx- The set of exponents {tx ,7xy} are not independent and scaling relations 



have the form [jT5|,|T4| 



= 1 + 

Ixy 

From (|3.71 ) one can find the simple expressions 



Ixy — 'Jxz'Jzy 



(3.7) 



(3.8) 



We have 10 unknown exponents altogether, namely and -jxy = 7^^, where x,y E {a, s, t, r}, 
but there exists only 6 linearly independent scaling relations ( p.7|) among them. Additional 
scaling relations can be obtained from the specific structure and evolution of an avalanche 
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and depend on the given model. The compactness of an avalanche cluster gives us 



2 16,14 



Thus, estimating only three critical exponents from the numerical data, we can calculate 
all the others using the scaling relations Eqs. (|3.71 ). Having calculated more than three 
exponents we are able to check these relations as well. The accurate determination of the 
r's exponents is a more difficult task than the 7's due to the strict dependence of the r's on 
the system size L. 

To reduce the fluctuations of the data, we integrated each distribution over bin lengths. 
The exponents j^y, x,y = {s, a, t}, are measured from the slopes of the straight parts of the 
corresponding graphs (Figs. ^). The obtained values are shown in Table |T[ 

Plotting integrated distributions P{t, L) -L^^ versus t-L'"^ on a double logarithmic scale, 
as is shown in Fig. ^ for different lattice sizes L, we obtained from finite-size scaling analysis 
that the best data collapse corresponds to (3t = 1.78 ± 0.05, = 1-36 ± 0.05 (Fig. ^). The 
scaling relation for the critical exponents (^) gives the value = 1.31 ± 0.05. 

Next, we use the measured values of r^, 7^4 and ^sa to estimate the whole set of exponents 
using the scaling relations, Eqs. (p.7[). These values are presented in Table fT[ 
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TABLES 



TABLE L The time averaged density pc of occupied sites for the HA model on the square 
lattice with two slightly different types of dynamics and on the Bethe lattice is compared with the 
value for the two-state Manna model. The uncertainty of the numerical data is about ±0.001. 









Model 




Density 


HA^ 


HA^ 


HA= 


Manna |§1 


Pc 


0.721 


0.755 


1 ^ 0.666 


0.683 



^N-E-S-W-N dynamics 
^N-S-W-E-N dynamics 
■^Bethe lattice 



TABLE n. The critical exponents for the 2d HA model evaluated in our work (first column) are 
compared with those for the BTW and two-state Manna models. The second column is the critical 
exponents for the BTW model obtained from numerical simulations, whereas in the third column 



we show exact values of the exponents for the BTW model based on the scaling relations (3.7), 
(|3.8| ) and jsr = Tr + 1 [16]. Comparison of the critical exponents of the HA and BTW models 
evaluated from numerical simulations shows that the HA model belongs to the universality class 
of the BTW model. The uncertainty of the numerical data for the HA model is about ±0.05. 

Model 



Exponent 


HA 


BTW 


BTW^ 


Manna 


Ts 


1.18^ 


1.20 


1 = 1.2 @ 


1.30 i 


Ta 


1.21^ 


1.22 [g 


! = 1-25 


1.37b 


Tt 


1.31 


1.32 


I = 1-4^ 


1.50 i 


Tr 


1.41^^ 


1.42^ 


i = 1-5^ 


1.75b 


Isa 


1.11 


1.06 [| 


f = 1.25^ 


1.23 i 


1st 


1.68 


1.64 
2.16^ 


2b 


1.67 i 


Isr 


2.23*' 


1 = 2.5b 


2.49b 


lat 


1.51 


1.52"^ 


1 = 1-6^ 


1.35 i 


lar 


2'= 


2= 


2" 


2.01b 


Itr 


1.33^ 


1.32 [| 


f = L25 


1.49 § 



^Exact result. 

bThe value of the exponent is obtained from the scaling relations ( |3.7D and ( |3.8| ). 
'^From the compactness of an avalanche cluster |16|. 



7 



The simulations for the HA model with N-S-W-E-N dynamics within a small uncertainty 
give the same values for the critical exponents. 



IV. THE HEIGHT- ARROW MODEL ON THE BETHE LATTICE. 

In this section, we present exact analytical calculations for the averaged density of occu- 
pied sites in the HA model on the Bethe lattice of coordination number q = 4. The Bethe 
lattice is defined through a Cayle tree well-known in graph theory which is a connected graph 
with no closed circuits of edges. Then, the Bethe lattice is an infinite Cayle tree homogenous 



in the sense that all except the outer vertices have the same coordination number q pO 



Following Dhar we approach the problem by dividing the allowed configurations of 
the HA model in the SOC state into two types: strongly allowed and weakly allowed and 
constructing the recurrent relations for the ratio of these configurations on the branches of 
the Bethe lattice. Using this ratio in the thermodynamic limit, we obtain the density of 
occupied sites in the HA model. 

First, let us briefly describe the procedure of construction of the Cayle tree. Like many 
tree-like structures, the Cayle tree of k generations of coordination number q can be con- 
structed by attaching q /cth-generation branches to a central site, as is shown in Fig. In 
turn, every A;th-generation branch is constructed by connecting q — 1 [k — l)th-generation 



branches to a new root and so on |20|. This property allows us to build recursion relations 
for the number of allowed configurations on the branches of the Cayle tree. 

The number of boundary sites of the Cayle tree is comparable with interior ones. Hence, 
the calculation of the bulk properties in the thermodynamic limit requires special care. Since 
we are interested in the solution on the Bethe lattice, we will take the result for the averaged 
density of occupied sites calculated at the central site of the Cayle tree as the value for the 
Bethe lattice. 

The definition of the HA model on this connected graph of coordination number g = 4 
remains unchanged. The only difference from the square lattice concerns the notation of 
directions of arrows and sequences of their turns. We will consider only sequential clockwise 
turns by the right angle and denote the directions of arrows at each site simply by {1, 2, 3, 4}. 

Let C be an allowed configuration on the A;th-generation branch Tk with root vertex a. 
Adding a vertex 6 to T^, one defines a subgraph T' = U b. If the sub configuration C on 
T' with Zh = and an arrow directed up or right (Fig. ^) becomes forbidden, C is called a 
weakly allowed (W) configuration, otherwise it is called a strongly allowed (S) one. 

Now consider Tj. with a root a that consists of three {k — l)th- generation branches 
T^^-i, ^fc^^ and Tj:^^ with roots oi, a2 and 03, respectively (Fig. ^). Let NY/{Tk,n, ]) and 
Ns{Tk,n,]) be the numbers of distinct W- and S'-type configurations on Tk with a given 
height Za = n and direction of the arrow at the root vertex a. 

Let us also introduce 

Nw{Tk)= E E Nw{Tk.n,r), (4.1) 

„={0,l}r={Ta} 

Ns{Tu)= E E Ns{n,n,r), (4.2) 

n={0,l}r={Ta} 
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where the first summation is over the values of the heights and the second one is over the 
directions of the arrow. As has been already mentioned, the arrow at each site can take only 
two directions. 

These numbers can be expressed in terms of the numbers of allowed subconfigurations 
on the three (A; — l)th-generation branches Tj^^i, T^li and T^li- 

Nwin) = iV<^)iVf ivf + N^^N^N^ + N^^^N^^N^^ + N^^^ N^^^ N^^ + 

iV«ivf ATf + iV«iVf Ar(^) + iV«iv£)ivf + Ar«ivg)iv(^) , (4.3) 

Ns{Tk) = 3N^s'^nPnP + 2N^^^nPnII) + 2N^^^ N^y^ + N^^^ N^^ N^^^ + 

2iV«ivf ATf + iV«Arf ivg) + Ar«iv£)ivf , (4.4) 

where A^W = N^{t1^1^), a = W,S and i = l, 2, 3. 
Let us define 

X = (4.5) 



If we consider graphs tI^\, Tjf]^ and T^^^i to be isomorphic, then A^(T^^i\) = N{tI^\^ 

.(3) 



A^(Tfc^_\) and from (|J) and (|J) one obtains the following recursion relation: 

X(Tfc) = ^(l + X(rfc_i)) . (4.6) 
With the initial condition X{To) = |, this equation has a simple solution 

X{n) = \-U-^'^'K (4.7) 

In the thermodynamic limit {k oo) the iterative sequence {X{Tk)} converges to a stable 
point X* = I that characterizes the ratio of the weakly allowed configurations to the strongly 
allowed ones in the SOC state. 

Consider now a randomly chosen site O deep inside the Cayle tree (Fig. |l^). The 
probability -P(l) of occupation of the site O is 

P{1) = §^, (4.8) 

where A^(l) is the number of allowed configurations with 2; = 1 at the site O and Atotai = 
A^(0) + A^(l) is the total number of allowed configurations on the Cayle tree. The numbers 
A^(0) and A^(l) can be expressed via the numbers of allowed configurations on the four 
neighbor /cth-generation branches T^^*\ i = 1,2,3,4 

4 

N{0) = 2[1 + 2X + X^]l[Ns{T^'^) , (4.9) 

i=l 

4 

iV(l) = 2[1 + AX + ^ 2X^] n NsiT^'^) . (4.10) 

i=l 
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For the sites far from the surface in the thermodynamic hmit {k — > oo) we have X = |. 
Thus, from( [4.9|) and ( 4.10 ) we obtain 



P(0) 



1 



P(l) 



2 



(4.11) 



3 



3 



The value for the concentration of occupied sites -P(l) is in good quahtative agreement with 
the numerical result obtained on the square lattice. 



We numerically studied the self-organizing height-arrow (HA) model on the square lattice 
and analytically on the Bethe lattice. The dynamics of the model drives it into the critical 
attractor with spatio-temporal complexity. The obtained distributions of various dynamic 
characteristics show an explicit power law behavior which indicates long-range correlations 
in the steady state of the system. To obtain the critical exponents of distributions of dynamic 
quantities of the model in the SOC state, we applied the finite-size scaling analysis. The 
values of these exponents are listed in Table |I| and compared with known exponents of the 
BTW model and two-state Manna model. Thus, we argue that the HA model belongs to 
the universality class of undirected models. 

Furthermore, we investigated the averaged density of occupied sites Pc in the SOC state 
of the HA model. It was also calculated exactly on the Bethe lattice of coordination number 
g = 4. The obtained results are presented in Table |. 
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FIGURES 



FIG. 1. A computer simulation of the HA model (N-E-S-W-N dynamics) on the square lattice 
of the linear size L = 600 with open boundary conditions. The steady state is reached by the 
model after about 100 000 avalanches. 

FIG. 2. The dependence of the time averaged density of occupied sites ((p{z = 1))l on the 
lattice size L. (a) The HA model with N-E-S-W-N dynamics and random initial directions of 
arrows, (b) The same model with N-S-W-E-N dynamics and arrows initially directed east or 
south. The numerical extrapolation for the infinity lattice size L — > co gives (a) Pc = 0.721 =b 0.001 
and (b) Pc = 0.755 ± 0.001. 

FIG. 3. A typical form of an avalanche cluster of the HA model. The lattice size is L = 200. 
The avalanche cluster has a layered structure. The number of topplings in each layer is indicated 
in gray scale. 

FIG. 4. Simulation results for distributions of avalanches by (a) mass (b) area and (c) duration 
of the HA model in the SOC state. The linear size of the lattice is L = 600. The number of 
avalanches for each distribution is 10^. 

FIG. 5. Double-logarithmic plot of the dependence of the stochastic variables {s,a,t} against 
each other for different lattice sizes. The distributions are integrated over bin lengths. 

FIG. 6. Double-logarithmic plot of the integrated distribution of avalanches P{t, L) versus du- 
ration t for five lattice sizes L = 100, 200, 500. Each distribution is averaged over 10^ avalanches. 

FIG. 7. Finite-size scaling plot for the integrated distributions P{t,L). The data for different 
L collapse onto a single curve for = 1.78 and i^t = 1.36 

FIG. 8. Construction of the Cayle tree with q = 4 and k = 3 generations by attaching g = 4 
/sth-generation branches to a central site. This procedure is explained in the text. 

FIG. 9. A fcth- generation branch and vertex b form a subgraph T'. The ovals denote the 
rest of the subbranches of T^. 

FIG. 10. A /cth-generation branch consists of three nearest {k — l)th-generation branches 
Ta,T<!',andT<?, 

FIG. 11. A site O with height Zo = n and a given direction of the arrow is located deep inside 
the lattice and surrounded by the four fcth-generation branches T^\t^\t^^ and . 
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FIG. 8. Construction of the Cayle tree with q = 4 and k = 3 generations by attaching q 
feth-generation branches to a central site. This procedure is explained in the text. 




FIG. 9. A A;th-generation branch and vertex b form a subgraph T'. The ovals denote the 
rest of the subbranches of T^. 




FIG. 10. A fcth-generation branch consists of three nearest {k — l)th-generation branches 

ii, 4^2, and Tt% 




FIG. 11. A site O with height Zq = n and a given direction of the arrow is located deep inside 
the lattice and surrounded by the four feth-generation branches T^^^ , xj^^^ , Tjf^ and T^^^ . 
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FIG. 1. A computer simulation of the HA model (N-E-S-W-N dynamics) 
on the square lattice of the linear size L=600 with open boundary conditions. 
The steady state is reached by the model after about 100000 avalanches. 




FIG. 2. The dependence of the time averaged density of occupied 
sites <p(z=l)> on the lattice size L. (a) The HA model with N-E-S-W-N 
dynamics and random initial directions of arrows, (b) The same model 
with N-S-W-E-N dynamics and arrows initially directed east or south. 
The numerical extrapolation for the infinity lattice size (L ^oo) gives 
(a) p^= 0.721 ±0.001 and (b) p^= 0.755 ±0.001. 




FIG. 2.(b) 




FIG. 3: A typical form of an avalanche cluster of the HA model. The lattice size is 
L = 200. The avalanche cluster has a layered structure. The number of topplings in each 
layer is indicated in gray scale. 




FIG. 4. Simulation results for distributions of avalanches by (a) mass, 
(b) area, and (c) duration of the HA model in the SOC state. The linear 
size of the lattice is L=600. The number of avalanches for each 
distribution is 10^. 




FIG. 4. (b). 




FIG. 4. (c). 




FIG. 5. Double-logarithmic plot of the dependence of the stochastic 
variables {s, a, t} against each other for different lattice sizes. 
The distributions are integrated over bin lengths. 




FIG.5.(b) 




FIG.5.(c) 
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FIG. 6. Double-logarithmic plot of the integrated distribution of 
avalanches P(t,L) versus duration t for the five lattice sizes 
L = 100, 200 500. Each distribution is averaged over 10' avalanches. 
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FIG. 7. Finite-size scaling plot for the integrated distributions P(t,L). 
The data for different L collapse onto a single curve for Pt = 1.78 and 
vt= 1.36. 



